#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@title: task2 

@author: dai
"""

import pylab as pl
import numpy as np
from radiosim import pyradio

def sine_est(sine_wave, fs=100, f_res = 0.1): 
    '''返回正弦波的幅度、频率、和相位，和投影向量'''  
    
    f_obs = np.arange(0, fs//2, f_res)
    N = sine_wave.shape[0]
    T = N/fs                  # total time
    prj_coef = []
    for fr in f_obs:
        x_ref = pyradio.sine_wave(ampl=1, fc=fr, phase=0, duration=T, fs=fs)[0]
        x_b = np.sum(sine_wave * np.conj(x_ref))/N  
        prj_coef.append(x_b)
    
    f_idx = np.argmax(np.abs(prj_coef))
    F_est = f_obs[f_idx]                       # 频率    
    A_est = prj_coef[f_idx] * np.sqrt(2)     # 幅度
    Phs_est = np.angle(A_est)                  # 相位  
    return A_est, F_est, Phs_est, prj_coef


#%% 观测信号
duration = 0.01
fs = 32000        # 采样频率
x_o, A1, fc1, phi1 = pyradio.rand_wave(duration=duration, fs=fs)


x_o = x_o.real  # 只取实部作为观测信号
# pyradio.plot_wave(x_o, fs);


#%% 参数估计与绘图
A, F, Phs, prj = sine_est(x_o, fs=fs, f_res = 10)   # 估计值
print(f'估计参数 幅度{np.abs(A):.4f}, 频率：{F:.4f}, 相位: {Phs:.4f}')

print(f'估计参数 幅度{A1:.4f}, 频率：{fc1:.4f}, 相位: {phi1:.4f}')

pl.plot(np.abs(prj))                                # 相关系数曲线
pl.grid()
pl.show()